其他
明码标价之任意科研图表绘制（以氨基酸的位点变异图为例）
前面我们开通了明码标价专栏：
都是公共数据的处理，其实也同步给了全部的代码，也算是一种粉丝福利吧！然后吸引到了各式各样的粉丝需求，仿佛都在这个月爆发了，不过我们人手有限，绝大部分需求只能说是牵线搭桥而已。比如有粉丝咨询科研绘图，我就会发在我们的兼职群，让有能力的人对接，比较简单就直接安排给学徒啦！大家可以在我们《生信技能树》公众号后台直接描述你的需求或者发邮件给到我哈（jmzeng1314@163.com）, 下面是一个例子，我安排给了学徒的！
最近曾老师让我画氨基酸的位点变异图，突然想到了maftools 中的🍭棒棒糖图：
lollipopPlot(maf = laml, gene = 'DNMT3A', AACol = 'Protein_Change', showMutationRate = TRUE)
但这样图还是获得不了位点的具体位置信息，图的可定制程度太低啦。
于是反手就是一个R 包：
trackViewer: a Bioconductor package for interactive and integrative visualization of multi-omics data. Nat Methods 16, 453–454 (2019) doi:10.1038/s41592-019-0430-y
trackViewer整合了Gviz
简单看了下：
画出来的图也听好看：
ps：感觉相关的知识点还是挺多的，太多懵懵懂懂，慢慢理解吧。
再前言
原谅我有些啰嗦。
这里主要是懵懵懂懂的尝试了这个包（或者说一个完整的新的类）之后的一些体悟。
如果只是针对trackViewer 画的棒棒糖图而言，其实是分成两个部分的：
一个是上方的mut/methy 部分，包含了样本的具体信息；一个是下方的block 也就是基因组部分，包含了参考的信息。
就像是R 的base 画图包一样，如果需要画某个基因的变异位点，完全可以自己创建Grange 对象，并调整对象后画图，当然这很麻烦；当然，对于变异数据，我们也可以导入vcf 数据，从中提取相关信息，进行绘图。
1. 自定义自己的棒棒糖图
1）绘制棒棒图
官方的描述是，棒棒糖图可以用来绘制methylation/variant/mutation 数据。
比如我们可以创建第一个自定义的棒棒糖图：
SNP <- c(10, 12, 23, 1402) # 指定SNP 坐标位置
sample.gr <- GRanges("chr1", IRanges(SNP, width=1, names=paste0("snp", SNP))) # 设置棒棒图位置
features <- GRanges("chr1", IRanges(c(1, 501, 1001), # 设置block起使位置
width=c(120, 400, 405), # 设置block 的长度
names=paste0("block", 1:3))) # 设置名字
lolliplot(sample.gr, features = features) # 绘图
我们还可以增加多一些SNP 标记：
SNP <- c(10, 100, 105, 108, 400, 410, 420, 600, 700, 805, 840, 1400, 1402)# 指定SNP 坐标位置
sample.gr <- GRanges("chr1", IRanges(SNP, width=1, names=paste0("snp", SNP))) # 设置棒棒图位置
features <- GRanges("chr1", IRanges(c(1, 501, 1001), # 设置block起使位置
width=c(120, 400, 405), # 设置block 的长度
names=paste0("block", 1:3))) # 设置名字
lolliplot(sample.gr, features = features) # 绘图
2）改棒棒糖颜色
利用RColorBrewer 获取颜色，为block 图像增加颜色：
# 2. painting ur plot
library(RColorBrewer)
features$fill <- brewer.pal(3, "Set1")
lolliplot(sample.gr, features)
除了改变block，还可以改变一下棒棒糖的“棒子”border 和棒棒糖的“糖”color，这里可以使用对象中默认的绘图颜色， 直接传入不同或相同的数字即可：
# candy
sample.gr$color <- sample.int(length(SNP), length(SNP))
# stick
sample.gr$border <- sample(c("gray80", "gray30"), length(SNP), replace=TRUE)
lolliplot(sample.gr, features)
我们还可以改变图像的透明度，这里的透明度应该是以0-1 为单位变化：
sample.gr2 <- sample.gr
sample.gr2$alpha <- sample(50:200, length(SNP), replace = TRUE)/200
lolliplot(sample.gr2, features)
3）添加索引及legend
sample.gr$label <- as.character(1:length(sample.gr))
sample.gr$label.col <- "black"
lolliplot(sample.gr, features)
我们可以添加一个legend：
# 为node 添加legend
my_color <- sample(c(brewer.pal(13, "Set3"), "#66C2A5"), 6)
sample.gr$color <- sample.int(6, 13, replace = T)
legend <- 1:6
names(legend) <- paste0("legend", letters[1:6])
lolliplot(sample.gr, features, legend=legend, xaxis=xaxis)
如果我们要对多个track或转录本绘棒棒糖图，最好使用列表定义legend。
4）改变block 高度
通过对feature 对象修改，我们可以获得坐标上的block，在创建时使用width 可以指定其长度，我们也可以后续通过设定height 属性改变它的高度：
features$height <- c(0.02, 0.05, 0.04)
lolliplot(sample.gr, features)
我们还可以用列表形式传递参数，并以数值, 单位的形式设置高度的单位：
features$height <- list(unit(1/16, "inches"),
unit(3, "mm"),
unit(12, "points"))
lolliplot(sample.gr, features)
5）改变转录本高度
通过属性score 进行调整：
# 改变node 高度/数目
sample.gr$score <- sample.int(5, length(sample.gr), replace = TRUE)
lolliplot(sample.gr, features)
6）绘制多个转录本
使用featureLayerID 属性控制多个转录本：
# 5. 绘制多个transcripts
features.mul <- rep(features, 2)
features.mul$fill <- brewer.pal(6, "Set1")
features.mul$featureLayerID <-
paste("tx", rep(1:2, each=length(features)), sep="_")
names(features.mul) <-
paste(features.mul$featureLayerID,
rep(1:length(features), 2), sep="_")
features.mul$height[4:6] <- list(unit(1/8, "inches"),
unit(0.5, "lines"),
unit(.2, "char"))
lolliplot(sample.gr, features.mul)
7）修改绘图函数参数
使用lolliplot 绘图，我们还可以直接对其参数进行修改：
lolliplot(
SNP.gr,
features = NULL,
ranges = NULL,
type = "circle",
newpage = TRUE,
ylab = TRUE,
ylab.gp = gpar(col = "black"),
yaxis = TRUE,
yaxis.gp = gpar(col = "black"),
xaxis = TRUE,
xaxis.gp = gpar(col = "black"),
legend = NULL,
cex = 1,
dashline.col = "gray80",
jitter = c("node", "label"),
rescale = FALSE,
...
)
Arguments
SNP.gr
A object of GRanges, GRangesList or a list of GRanges. All the width of GRanges must be 1.
features
A object of GRanges, GRangesList or a list of GRanges. The metadata 'featureLayerID' are used for drawing features in different layers. See details in vignette.
ranges
A object of GRanges or GRangesList.
type
character. Could be circle, pie, pin, pie.stack or flag.
newpage
Plot in the new page or not.
ylab
Plot ylab or not. If it is a character vector, the vector will be used as ylab.
ylab.gp, xaxis.gp, yaxis.gp
An object of class gpar for ylab, xaxis or yaxis.
yaxis
Plot yaxis or not.
xaxis
Plot xaxis or not. If it is a numeric vector with length greater than 1, the vector will be used as the points at which tick-marks are to be drawn. And the names of the vector will be used to as labels to be placed at the tick points if it has names.
legend
If it is a list with named color vectors, a legend will be added.
cex
cex will control the size of circle.
dashline.col
color for the dashed line.
jitter
jitter the position of nodes or labels.
rescale
logical(1) or a dataframe with rescale from and to. Recalse the x-axis or not. if dataframe is used, colnames must be from.start, from.end, to.start, to.end.
比如可以去掉坐标轴：
lolliplot(sample.gr, features, yaxis = F, xaxis = F)
比如将circle 风格替换为pin ：
lolliplot(sample.gr, features, type="pin")
比如flag 风格：
# 国旗型
sample.gr.flag <- sample.gr
sample.gr.flag$label <- names(sample.gr) ## move the names to metadata:label
names(sample.gr.flag) <- NULL
lolliplot(sample.gr.flag, features, type="flag", ranges=GRanges("chr1", IRanges(0, 1600)))
饼图，比如也可以按照饼图的比例表示不同类型突变数目的比值情况：
# 饼图
sample.gr$score <- NULL
sample.gr$label <- NULL
# 清空所有其他无关数值
sample.gr$color <- rep(list(c("#87CEFA", '#98CE31')), length(SNP)) # 这里需要使用列表，对应饼图各部分颜色
sample.gr$border <- "gray30"
sample.gr$value1 <- sample.int(5, length(sample.gr), replace = TRUE) # 饼图数值1 所占比重
sample.gr$value2 <- 10 - sample.gr$value1 # 数值2 所占比重
lolliplot(sample.gr, features, type="pie", xaxis=xaxis)
我们还可以结合先前的内容，加上legend：
8）调整x与y 轴
# 6. 调整x 与y 轴
xaxis <- c(1, 200, 400, 701, 1000, 1200, 1402) ## define the position
names(xaxis) <- xaxis
# names(xaxis)[4] <- "center"
lolliplot(sample.gr, features, xaxis=xaxis)
通过命令向量，我们可以修改legand 上的内容，比如将中间命名为center：
# 6. 调整x 与y 轴
xaxis <- c(1, 200, 400, 701, 1000, 1200, 1402) ## define the position
names(xaxis) <- xaxis
names(xaxis)[4] <- "center"
lolliplot(sample.gr, features, xaxis=xaxis)
y 轴同理。
9）修改label 与title 文本
该绘图函数推测是基于基础包开发，因此很多图形参数可以使用基础包的参数或函数：
# 7. 修改title 与label text
sample.gr$score <- sample.int(5, length(sample.gr), replace = TRUE)
names(xaxis)[4] <- "center"
sample.gr$label.parameter.rot <- 60
# ylab
lolliplot(sample.gr, features, legend=legend, ylab="mut counts")
# xlab
grid.text("location", x=.5, y=.01, just="bottom")
# title
grid.text("SUPER: [MUTATION RATE:66.666%]", x=.5, y=.98, just="top",
gp=gpar(cex=1.5, fontface="bold"))
2. 载入VCF 数据
2.0 载入相关包
library(VariantAnnotation)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)
library(trackViewer)
2.1 选择特定基因
这里就选择文件中出现的第一个基因OR4F5 尝试一下https://www.genecards.org/cgi-bin/carddisp.pl?gene=OR4F5 ：
首先创建一段基因range，这里我选择一号染色体的69200, 70000 （基本涵盖了上面的起始与终止位点）：
my_gr <- GRanges("chr1", IRanges(69200, 70000, names="OR4F5"))
2.2 读取vcf 文件
接着需要读取vcf文件，但我发现如果直接读取会报错：
> dir(path = annovar_dir)[2]
[1] "annovarhg38_multianno.vcf"
> vcf <- readVcf(annovar_d_file, "hg38", param=my_gr) # param 指定读取的片段
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘scanVcf’ for signature ‘"character", "GRanges"’
思考了一下，发现需要对文件索引创建一个对象，而注释后的vcf 文件并没有索引，另外，范例使用的是gz 格式数据，为了加速函数读取，也对annovar 注释结果gz压缩一下（这里在linux 完成）：
# gz 压缩
$ bgzip -c annovarhg38_multianno.vcf > annovarhg38_multianno.vcf.gz
# 使用-t 创建tbi 索引
bcftools index -t annovarhg38_multianno.vcf.gz
我们可以通过VariantAnnotation 包中的readVcf 函数读取指定片段的vcf 信息：
> annovar_d_file
# 读取vcf 索引文件
tab <- TabixFile(annovar_d_file)
# 读取vcf 文件
vcf <- readVcf(annovar_d_file, "hg38", param=my_gr) # param 指定读取的片段
2.3 提取数据
首先是变异数据位点信息以及具体的信息列合并：
mutation.gr <- rowRanges(vcf)
# mcol 可以合并S4 对象中的表格与一般df
# 获得vcf 文件的信息列,合并位点记录信息与信息列
mcols(mutation.gr) <- cbind(mcols(mutation.gr),
VariantAnnotation::info(vcf))
# 获得AF百分比%
mutation.gr$score <- mutation.gr$AF*100
因为数据库使用的是UCSC，因此还需要将其转换一下：
# 修改下基因类型，注释也会使用NCBI，切换为UCSC
# seqlevelsStyle(mutation.gr) <- "UCSC"
seqlevelsStyle(gr) <- "UCSC"
trs <- geneModelFromTxdb(TxDb.Hsapiens.UCSC.hg38.knownGene,
org.Hs.eg.db,
gr=my_gr)
2.4 绘图
# 准备绘图
library(RColorBrewer)
## 棒棒糖
## 先修改mutation.gr 名称
# names(mutation.gr)
names(mutation.gr) <- sub("chr1:", '', names(mutation.gr))
names(mutation.gr)
## 设置颜色
mutation.gr$border <- "gray30"
## 按突变类型修改颜色
#mutation.gr$ExonicFunc.refGene
mutation.gr$color <-
ifelse(grepl("nonsynonymous_SNV", mutation.gr$ExonicFunc.knownGene), "lightcyan", ifelse(grepl("synonymous_SNV", mutation.gr$ExonicFunc.knownGene), "lavender", "#FFFFB3"))
## 添加index
mutation.gr$label <- as.character(1:length(mutation.gr))
## 添加legend
legends <- list(labels=c("nonsynonymous", "synonymous", "stopgain"), fill=c("lightcyan", 'lavender', '#FFFFB3'))
## block
features <- trs[[1]]$dat
features$fill <- brewer.pal(3, "Set1")
names(features) <- features$feature
## 画图
## y轴加标记
lolliplot(mutation.gr, features = features, legend = legends,
ylab="Mut Freq%")
## 标题
grid.text("OR4F5", x=.5, y=.90, just="top",
gp=gpar(cex=1.5, fontface="bold"))
## x轴标记
grid.text("chr1", x=0.9, y=.07, just="bottom", gp=gpar(cex=1.5, fontface="bold"))
其他知识
Grange/IRanges 对象
Grange/IRanges 对象用来储存染色体的相关信息，主要有三个属性（name、seqnames染色体名、ranges长度、strand 正or负链），以及一些其他的属性比如绘图属性（fill, height）：
> features.mul
GRanges object with 6 ranges and 2 metadata columns:
seqnames ranges strand | fill height
<Rle> <IRanges> <Rle> | <character> <numeric>
block1 chr1 1-120 * | #E41A1C 0.02
block2 chr1 501-900 * | #377EB8 0.05
block3 chr1 1001-1405 * | #4DAF4A 0.04
直接查看该对象，也可以看到不同类型属性的分割线：
创建时我们可以进行指定：
features <- GRanges("chr1", IRanges(c(1, 501, 1001), # 设置block起使位置
width=c(120, 400, 405), # 设置block 的长度
names=paste0("block", 1:3))) # 设置名字
默认下strand 为+/* 我们可以用strand 函数定义修改：
> strand(features) <- "-"
> features
GRanges object with 3 ranges and 2 metadata columns:
seqnames ranges strand | fill height
<Rle> <IRanges> <Rle> | <character> <numeric>
block1 chr1 1-120 - | #E41A1C 0.02
block2 chr1 501-900 - | #377EB8 0.05
block3 chr1 1001-1405 - | #4DAF4A 0.04
start(), end(), width() 函数可以获得Grange/IRanges 对象位点的起使、终止、全部的长度信息。
对于分段的小的基因片段：
基因 TYMP
trs[[1]]$dat
GRanges object with 11 ranges and 7 metadata columns:
seqnames ranges strand | feature id exon
<Rle> <IRanges> <Rle> | <character> <character> <character>
[1] chr22 50964181-50964198 - | utr3 unknown uc003bmb.5_9
[2] chr22 50964199-50964347 - | CDS unknown uc003bmb.5_9
[3] chr22 50964430-50964570 - | CDS unknown uc003bmb.5_8
[4] chr22 50964675-50964905 - | CDS unknown uc003bmb.5_7
[5] chr22 50965005-50965167 - | CDS unknown uc003bmb.5_6
[6] chr22 50965594-50965712 - | CDS unknown uc003bmb.5_5
[7] chr22 50966017-50966146 - | CDS unknown uc003bmb.5_4
[8] chr22 50966941-50967039 - | CDS unknown uc003bmb.5_3
[9] chr22 50967565-50967767 - | CDS unknown uc003bmb.5_2
[10] chr22 50967925-50968138 - | CDS unknown uc003bmb.5_1
[11] chr22 50968139-50968514 - | utr5 unknown uc003bmb.5_1
transcript gene symbol density
<character> <character> <character> <numeric>
[1] uc003bmb.5 1890 TYMP 1
[2] uc003bmb.5 1890 TYMP 1
[3] uc003bmb.5 1890 TYMP 1
[4] uc003bmb.5 1890 TYMP 1
[5] uc003bmb.5 1890 TYMP 1
[6] uc003bmb.5 1890 TYMP 1
[7] uc003bmb.5 1890 TYMP 1
[8] uc003bmb.5 1890 TYMP 1
[9] uc003bmb.5 1890 TYMP 1
[10] uc003bmb.5 1890 TYMP 1
[11] uc003bmb.5 1890 TYMP 1
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
基因 ODF3B：
trs[[5]]$dat
GRanges object with 8 ranges and 7 metadata columns:
seqnames ranges strand | feature id exon
<Rle> <IRanges> <Rle> | <character> <character> <character>
[1] chr22 50968838-50968870 - | utr3 unknown uc003bmg.2_6
[2] chr22 50968871-50968989 - | CDS unknown uc003bmg.2_6
[3] chr22 50969141-50969281 - | CDS unknown uc003bmg.2_5
[4] chr22 50969366-50969487 - | CDS unknown uc003bmg.2_4
[5] chr22 50969620-50969724 - | CDS unknown uc003bmg.2_3
[6] chr22 50969999-50970131 - | CDS unknown uc003bmg.2_2
[7] chr22 50970398-50970506 - | CDS unknown uc003bmg.2_1
[8] chr22 50970507-50970540 - | utr5 unknown uc003bmg.2_1
transcript gene symbol density
<character> <character> <character> <numeric>
[1] uc003bmg.2 440836 ODF3B 1
[2] uc003bmg.2 440836 ODF3B 1
[3] uc003bmg.2 440836 ODF3B 1
[4] uc003bmg.2 440836 ODF3B 1
[5] uc003bmg.2 440836 ODF3B 1
[6] uc003bmg.2 440836 ODF3B 1
[7] uc003bmg.2 440836 ODF3B 1
[8] uc003bmg.2 440836 ODF3B 1
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
我们可以将这两个GRanges 对象整合起来，使用range 获得整段长度：
> range(trs[[1]]$dat)
GRanges object with 1 range and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr22 50964181-50968514 -
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
接着可以用向量连接它们：
features <- c(range(trs[[1]]$dat), range(trs[[5]]$dat))
> features
GRanges object with 2 ranges and 1 metadata column:
seqnames ranges strand | fill
<Rle> <IRanges> <Rle> | <character>
TYMP chr22 50964181-50968514 - | #E41A1C
ODF3B chr22 50968838-50970540 - | #377EB8
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
新的feature 就包含了两段基因的信息。
载入基因组信息
readVcf 函数可以用来读取vcf文件，可以通过设定参数param 指定一个Grange 对象，来对染色体指定区域的突变信息进行读取：
library(VariantAnnotation)
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
gr <- GRanges("22", IRanges(50968014, 50970514, names="TYMP"))
# 如果有索引，可以读取vcf 索引文件，加速读取
tab <- TabixFile(fl)
# 读取vcf 文件
vcf <- readVcf(fl, "hg19", param=gr) # param 指定读取的片段
网页工具完成不了如此个性化的图表哦
如果你没有时间从头开始学编程，也可以委托给我们这样的专业的团队，付费拿到同样的数据分析， 绝大部分科研图表绘制，都是800元哦，可以拿到全部的数据和代码哦！
需要自己读文献筛选合适的数据集，解读数据要求 提供半个小时左右的一对一讲解数据和代码
如果需要委托，直接在我们《生信技能树》公众号留言即可，我们会安排合适的生信工程师对接具体的项目。